{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Baseline examples\n",
    "\n",
    "Charles Le Losq\n",
    "\n",
    "August 2017\n",
    "\n",
    "Examples of using the baseline function with various algorithms (splines, polynomials, ALS, arPLS...)\n",
    "\n",
    "## Documentation\n",
    "\n",
    "See http://charlesll.github.io/Spectra.jl/stable/ for information about spectra.\n",
    "\n",
    "## References\n",
    "\n",
    "a good read about the ALS and arPLS algorithms is \n",
    "\n",
    "    Baek, S.-J., A. Park, Y.-J. Ahn, and J. Choo (2015), Baseline correction using asymmetrically reweighted penalized least squares smoothing, Analyst, 140(1), 250–257, doi:10.1039/C4AN01061B.\n",
    "\n",
    "## Importing the relevant libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The Julia library to treat Raman spectroscopy data\n",
    "using Spectra\n",
    "\n",
    "# to plot stuffs\n",
    "using PyPlot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a fake signal to know the ground truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAG9CAYAAADuqYmYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlYVGX/P/D3DJvIKioiAm65oCGY5p5pUu5WaqWR+ZjptxQlLXMhs6xHyzRLsawnyyzNpcXUytJc0FxCccMFN3JD3JBd1rl/f/A7pznMDAwwywHer+viauacMzP3kMB7PvemEUIIEBEREZFMa+8GEBEREakNAxIRERFRCQxIRERERCUwIBERERGVwIBEREREVAIDEhEREVEJDEhEREREJTAgEREREZXAgERERERUAgMSERERUQkMSEREREQlONq7AVWFTqdDcnIyPDw8oNFo7N0cIiIiMoMQApmZmfD394dWa35diAHJTMnJyQgMDLR3M4iIiKgCrly5goCAALOvZ0Ayk4eHB4Dib7Cnp6edW0NERETmyMjIQGBgoPx33FwMSGaSutU8PT0ZkIiIiKqY8g6P4SBtIiIiohIYkIiIiIhKYEAiIiIiKoFjkIiIiKxIp9MhPz/f3s2otpycnODg4GDx52VAIiIispL8/HwkJSVBp9PZuynVmre3N/z8/Cy6TiEDEhERkRUIIXD9+nU4ODggMDCwXIsUknmEEMjJycHNmzcBAA0bNrTYczMgERERWUFhYSFycnLg7++P2rVr27s51ZarqysA4ObNm/D19bVYdxvjLBERkRUUFRUBAJydne3ckupPCqAFBQUWe04GJCIiIivi/p3WZ43vMQMSERERUQkMSERERGSW//znP3jiiSds/rpvvfUWwsLCbPqaHKRNREREZvn4448hhLB3M2yCAYmIVC8nJweurq4cy0FkZ15eXvZugs2wi42IVO3UqVPw9PREVFSUvZtCVGN8//33CAkJgaurK+rWrYvw8HBkZ2cbdLFlZmYiIiICbm5uaNiwIRYvXoxevXrhlVdeka9p0qQJ5s2bhxdeeAEeHh4ICgrC559/rni96dOno2XLlqhduzaaNWuG2bNnW3RGWkWwgkREqhYdHY2ioiIsXboUS5YssXdziCpMWtTQHmrXrm12Bfb69esYOXIkFixYgCeffBKZmZnYs2eP0a61qVOn4q+//sKmTZvQoEEDvPnmm4iPjzcYL7Ro0SK88847mDVrFr7//nu8/PLLePjhh9GqVSsAgIeHB1auXAl/f3+cOHEC48aNg4eHB15//fXKv/kKYkAiIlW7deuWvZtAZBE5OTlwd3e3y2tnZWXBzc3NrGuvX7+OwsJCDB06FI0bNwYAhISEGFyXmZmJr7/+GmvWrEGfPn0AAF999RX8/f0Nrh0wYAAmTJgAoLhatHjxYuzcuVMOSG+88YZ8bZMmTfDaa69h7dq1DEhERKZIWwgQkW2EhoaiT58+CAkJQd++ffHYY49h+PDhqFOnjuK6ixcvoqCgAJ06dZKPeXl5yaFHX7t27eTbGo0Gfn5+ip/tdevWYcmSJbhw4QKysrJQWFgIT09PK7w78zEgEZGqMSBRdVG7dm1kZWXZ7bXN5eDggG3btmHfvn34448/sHTpUkRHR+PgwYMVfn0nJyfFfY1GI2/gu3//fkRERODtt99G37594eXlhbVr12LRokUVfj1LYEAiIlVLT0+3dxOILEKj0ZjdzWVvGo0G3bt3R/fu3fHmm2+icePG+OmnnxTXNGvWDE5OToiLi0NQUBCA4p/Xs2fPomfPnma/1r59+9C4cWNER0fLxy5dumSZN1IJDEhEREQkO3jwIP7880889thj8PX1xcGDB3Hr1i0EBwfj+PHj8nUeHh4YPXo0pk2bBh8fH/j6+mLOnDnQarXlWpKjRYsWuHz5MtauXYsHH3wQv/zyi0EYswdO8yci1bJXdwRRTebp6YnY2FgMGDAALVu2xBtvvIFFixahf//+Btd++OGH6Nq1KwYNGoTw8HB0794dwcHBqFWrltmvN2TIEEyZMgWRkZEICwvDvn37MHv2bEu+pQrRiJqyJGYlZWRkwMvLC+np6XYfOEZUU5w5cwbBwcEAindEz8vLs3OLiMyXm5uLpKQkNG3atFyBoSrLzs5Go0aNsGjRIowdO9Zmr1va97qif7/ZxUZEqnXlyhX5dn5+PoQQXE2bSEWOHDmCM2fOoFOnTkhPT8fcuXMBAI8//ridW1Z5DEhEpFrHjh1T3C8oKICzs7OdWkNExixcuBCJiYlwdnZGhw4dsGfPHtSrV8/ezao0BiQiUq2//vpLcT8vL48BiUhF2rdvj8OHD9u7GVbBQdpEpEpCCOzbt09xjGOQiMhWGJCISJUuXLiAmzdvKipGDEhUFXEulPVZ43vMgEREqiSttxIaGiovrseARFWJg4MDgOIJBmRd0ibAJVfsrgyOQSIiVcrMzAQA1K1bF87OzsjOzmZAoirF0dERtWvXxq1bt+Dk5AStljUJSxNCICcnBzdv3oS3t7ccSi2BAYmIVElaJNLd3R0uLi4A+EmcqhaNRoOGDRsiKSlJFVtnVGfe3t7w8/Oz6HMyIBGRKhkLSKwgUVXj7OyMFi1aMNxbkZOTk0UrRxIGJCJSJQYkqi60Wm2NWUm7OlF9h2hsbCwGDx4Mf39/aDQabNy4UXE+KysLkZGRCAgIgKurK9q0aYPly5crrsnNzcXEiRNRt25duLu7Y9iwYbhx44Yt3wYRlZMUkNzc3BiQiMjmVB+QsrOzERoaimXLlhk9P3XqVGzduhXffvstTp8+jVdeeQWRkZHYtGmTfM2UKVOwefNmbNiwAbt370ZycjKGDh1qq7dARBWgX0GSpvozIBGRrai+i61///5GdxCW7Nu3D6NHj0avXr0AAOPHj8dnn32Gv//+G0OGDEF6ejpWrFiBNWvW4JFHHgEAfPXVVwgODsaBAwfQpUsXW7wNIiqn7OxsAOxiIyL7UH0FqSzdunXDpk2bcO3aNQghsHPnTpw9exaPPfYYAODw4cMoKChAeHi4/JjWrVsjKCgI+/fvN/m8eXl5yMjIUHwRke1wFhsR2VOVD0hLly5FmzZtEBAQAGdnZ/Tr1w/Lli1Dz549AQApKSlwdnaGt7e34nENGjRASkqKyeedP38+vLy85K/AwECrvg8iUuIgbSKyp2oRkA4cOIBNmzbh8OHDWLRoESZOnIjt27dX6nlnzpyJ9PR0+evKlSsWajERmYMBiYjsSfVjkEpz7949zJo1Cz/99BMGDhwIAGjXrh2OHj2KhQsXIjw8HH5+fsjPz0daWpqiinTjxo1SF5VycXGRfykTke0xIBGRPVXpClJBQQEKCgoMlm93cHCATqcDAHTo0AFOTk74888/5fOJiYm4fPkyunbtatP2EpH5OIuNiOxJ9RWkrKwsnD9/Xr6flJSEo0ePwsfHB0FBQXj44Ycxbdo0uLq6onHjxti9ezdWrVqFDz/8EADg5eWFsWPHYurUqfDx8YGnpycmTZqErl27cgYbkYoZWweJg7SJyFZUH5AOHTqE3r17y/enTp0KABg9ejRWrlyJtWvXYubMmYiIiEBqaioaN26M//73v3jppZfkxyxevBharRbDhg1DXl4e+vbti08++cTm74WIzMcuNiKyJ9UHpF69ekEIYfK8n58fvvrqq1Kfo1atWli2bJnJxSaJSF3y8/NRUFAAgAGJiOyjSo9BIqLqSVokEuBWI0RkHwxIRKQ6Uveas7Oz/AUwIBGR7TAgEZHq6I8/AsAKEhHZHAMSEamOqYDEWWxEZCsMSESkOqwgEZG9MSARkeror4EEMCARke0xIBGR6rCCRET2xoBERKojTfOXAhJnsRGRrTEgEZHqlKwg1apVCwADEhHZDgMSEalOyYBUu3ZtAEBOTo7d2kRENQsDEhGpjqmApL/CNhGRNTEgEZHqlAxI0mw2BiQishUGJCJSHQYkIrI3BiQiUh1TAYljkIjIVhiQiEh1Si4UKf23sLCQ240QkU0wIBGR6pRcB0kapK1/jojImhiQiEh1SnaxOTs7w9HREQADEhHZBgMSEalOyYAEcKA2EdkWAxIRqU5pAYkDtYnIFhiQiEh1WEEiIntjQCIiVdHpdAaDtAGupk1EtsWARESqcu/ePQghALCCRET2w4BERKoida8BgKurq3ybAYmIbIkBiYhURQpAbm5u0Gr//RXFQdpEZEsMSESkKsYGaAMcg0REtsWARESqYiogsYuNiGyJAYmIVIUBiYjUgAGJiFSl5Ea1EgYkIrIlBiQiUpV79+4BUG5QC3CQNhHZFgMSEalKbm4uAKBWrVqK4xykTUS2xIBERKpiKiCxi42IbIkBiYhUhQGJiNRA9QEpNjYWgwcPhr+/PzQaDTZu3GhwzenTpzFkyBB4eXnBzc0NDz74IC5fviyfz83NxcSJE1G3bl24u7tj2LBhuHHjhi3fBhGZSRqDpL+KNvBvFxvHIBGRLag+IGVnZyM0NBTLli0zev7ChQvo0aMHWrdujV27duH48eOYPXu24tPnlClTsHnzZmzYsAG7d+9GcnIyhg4daqu3QETlYKqCJAUm6TwRkTU52rsBZenfvz/69+9v8nx0dDQGDBiABQsWyMeaN28u305PT8eKFSuwZs0aPPLIIwCAr776CsHBwThw4AC6dOlivcYTUbmZCkjSfQYkIrIF1VeQSqPT6fDLL7+gZcuW6Nu3L3x9fdG5c2dFN9zhw4dRUFCA8PBw+Vjr1q0RFBSE/fv3m3zuvLw8ZGRkKL6IyPpMdbExIBGRLVXpgHTz5k1kZWXhvffeQ79+/fDHH3/gySefxNChQ7F7924AQEpKCpydneHt7a14bIMGDZCSkmLyuefPnw8vLy/5KzAw0KrvhYiKsYuNiNSgSgcknU4HAHj88ccxZcoUhIWFYcaMGRg0aBCWL19eqeeeOXMm0tPT5a8rV65YoslEVIayutikChMRkTWpfgxSaerVqwdHR0e0adNGcTw4OBh79+4FAPj5+SE/Px9paWmKKtKNGzfg5+dn8rldXFzg4uJinYYTkUnsYiMiNajSFSRnZ2c8+OCDSExMVBw/e/YsGjduDADo0KEDnJyc8Oeff8rnExMTcfnyZXTt2tWm7SWispVVQcrPz5erx0RE1qL6ClJWVhbOnz8v309KSsLRo0fh4+ODoKAgTJs2Dc888wx69uyJ3r17Y+vWrdi8eTN27doFAPDy8sLYsWMxdepU+Pj4wNPTE5MmTULXrl05g41IhcoKSEDxJIqSFSYiIktSfUA6dOgQevfuLd+fOnUqAGD06NFYuXIlnnzySSxfvhzz58/H5MmT0apVK/zwww/o0aOH/JjFixdDq9Vi2LBhyMvLQ9++ffHJJ5/Y/L0QUdmkgFQyAOnfz83NZUAiIqvSCCGEvRtRFWRkZMDLywvp6enw9PS0d3OIqq327dvj6NGj2Lp1K/r27as45+joiKKiIly7dg3+/v52aiERVSUV/ftdpccgEVH1Y6qLTf8YB2oTkbUxIBGRqjAgEZEaMCARkaqYmuYPMCARke0wIBGRqpRWQeJq2kRkKwxIRKQq5nSxcTVtIrI2BiQiUg2dToe8vDwA7GIjIvtiQCIi1ZDCEcBB2kRkXwxIRKQa+sGHAYmI7IkBiYhUQwo+Dg4OcHJyMjgvdbtxDBIRWRsDEhGphhR8jFWP9I+zgkRE1saARESqUdoMNv3jDEhEZG0MSESkGgxIRKQWDEhEpBqlraINMCARke0wIBGRapRVQeIgbSKyFQYkIlINdrERkVowIBGRajAgEZFaMCARkWowIBGRWjAgEZFqSFuNuLi4GD3PgEREtsKARESqwUHaRKQWDEhEpBqsIBGRWjAgEZFqSMGnrIDEChIRWRsDEhGphlRBMtXF5uXlBQBIT0+3WZuIqGZiQCIi1Siri83b2xsAkJaWZrM2EVHNxIBERKpR1iBtBiQishUGJCJSDXMrSPfu3ZOvJSKyBgYkIlKNsgKSp6enfJvjkIjImhiQiEg1yupic3BwkEMSu9mIyJoYkIhINcqqIAEch0REtsGARESqUVYFCQDq1KkDgAGJiKyLAYmIVIMVJCJSCwYkIlINBiQiUgsGJCJSDXO62BiQiMgWVB+QYmNjMXjwYPj7+0Oj0WDjxo0mr33ppZeg0Wjw0UcfKY6npqYiIiICnp6e8Pb2xtixY5GVlWXtphNROZWngnT37l2btImIaibVB6Ts7GyEhoZi2bJlpV73008/4cCBA/D39zc4FxERgZMnT2Lbtm3YsmULYmNjMX78eGs1mYgqqKzNagFWkIjINhzt3YCy9O/fH/379y/1mmvXrmHSpEn4/fffMXDgQMW506dPY+vWrYiLi0PHjh0BAEuXLsWAAQOwcOFCo4GKiOyjrM1qAQYkIrIN1VeQyqLT6TBq1ChMmzYNbdu2NTi/f/9+eHt7y+EIAMLDw6HVanHw4EFbNpWIysBB2kSkFqqvIJXl/fffh6OjIyZPnmz0fEpKCnx9fRXHHB0d4ePjg5SUFJPPm5eXp9jrKSMjwzINJiKTzBmk7eXlBYBbjRCRdVXpCtLhw4fx8ccfY+XKldBoNBZ97vnz58PLy0v+CgwMtOjzE5EhcypIUniSwhQRkTVU6YC0Z88e3Lx5E0FBQXB0dISjoyMuXbqEV199FU2aNAEA+Pn54ebNm4rHFRYWIjU1FX5+fiafe+bMmUhPT5e/rly5Ys23QlTjCSHMCkjSufz8fJu0i4hqpirdxTZq1CiEh4crjvXt2xejRo3CmDFjAABdu3ZFWloaDh8+jA4dOgAAduzYAZ1Oh86dO5t8bhcXl1J/SRORZRUWFkKn0wEovYtN+rnU7wInIrI01QekrKwsnD9/Xr6flJSEo0ePwsfHB0FBQahbt67ieicnJ/j5+aFVq1YAgODgYPTr1w/jxo3D8uXLUVBQgMjISIwYMYIz2IhURD/wlPbhxNnZ2eB6IiJLU30X26FDh9C+fXu0b98eADB16lS0b98eb775ptnPsXr1arRu3Rp9+vTBgAED0KNHD3z++efWajIRVYD+mCJ2sRGRvam+gtSrVy8IIcy+/p9//jE45uPjgzVr1liwVURkaVJFyNHREQ4ODiavYxcbEdmC6itIRFQzmDNAG/i3i40VJCKyJgYkIlIFc9ZAAlhBIiLbYEAiIlUobwWpqKgIRUVFVm8XEdVMDEhEpArlrSAB7GYjIuthQCIiVTC3gqR/nt1sRGQtDEhEpArmBiQnJyeDxxARWRoDEhGpgrldbBqNhjPZiMjqGJCISBWkgGTOFj+cyUZE1saARESqcO/ePQCAq6trmdeygkRE1saARESqIFWQzAlIrCARkbUxIBGRKpSngsSARETWxoBERKrALjYiUhMGJCJSBSkglTWLDWAFiYisjwGJiFShIhUkBiQishYGJCJShYqMQWIXGxFZCwMSEakCZ7ERkZowIBGRKrCLjYjUhAGJiFSBXWxEpCYMSESkCpzFRkRqwoBERKrAdZCISE0YkIhIFbiSNlHFCCHs3YRqiQGJiFSBs9iIym/06NFo2bIlcnJy7N2UaocBiYhUgV1sROX3yy+/4Pz58zh37py9m1LtMCARkSqwi42o/AoKChT/JcthQCIiVSjPLDaug0RUjAHJehiQiEgVuA4SUfkxIFkPAxIR2Z0Qgl1sROUkhEBhYSEABiRrYEAiIrvTrwRxqxEi8xQVFcm3GZAsjwGJiOxOqh4B7GIjMpd+KGJAsjwGJCKyOykgabVaODk5lXk9u9iIGJCsjQGJiOxOfwabRqMp83qug0QEefwRwIBkDQxIRGR35RmgDbCCRASwgmRtDEhEZHflDUjSWkn6Y5eIahr9UMRqquWpPiDFxsZi8ODB8Pf3h0ajwcaNG+VzBQUFmD59OkJCQuDm5gZ/f388//zzSE5OVjxHamoqIiIi4OnpCW9vb4wdOxZZWVm2fitEZEJ59mEDAA8PDwBAZmam1dpEpHasIFmX6gNSdnY2QkNDsWzZMoNzOTk5iI+Px+zZsxEfH48ff/wRiYmJGDJkiOK6iIgInDx5Etu2bcOWLVsQGxuL8ePH2+otEFEZyltBYkAiYkCyNkd7N6As/fv3R//+/Y2e8/LywrZt2xTHYmJi0KlTJ1y+fBlBQUE4ffo0tm7diri4OHTs2BEAsHTpUgwYMAALFy6Ev7+/1d8DEZWOAYmo/DhI27pUX0Eqr/T0dGg0Gnh7ewMA9u/fD29vbzkcAUB4eDi0Wi0OHjxo8nny8vKQkZGh+CIi68jJyQFQ/oCUlZUFnU5ntXYRqRkrSNZVrQJSbm4upk+fjpEjR8LT0xMAkJKSAl9fX8V1jo6O8PHxQUpKisnnmj9/Pry8vOSvwMBAq7adqCbLzs4GALi7u5t1vRSQ9B9LVNMwIFlXtQlIBQUFePrppyGEwKefflrp55s5cybS09PlrytXrliglURkjDRpws3NzazrXV1dodUW//piNxvVVAxI1qX6MUjmkMLRpUuXsGPHDrl6BAB+fn64efOm4vrCwkKkpqbCz8/P5HO6uLjIa60QkXWVt4Kk0Wjg4eGB9PR0BiSqsTgGybqqfAVJCkfnzp3D9u3bUbduXcX5rl27Ii0tDYcPH5aP7dixAzqdDp07d7Z1c4nICKmCZG5AAjhQm4gVJOtSfQUpKysL58+fl+8nJSXh6NGj8PHxQcOGDTF8+HDEx8djy5YtKCoqkscV+fj4wNnZGcHBwejXrx/GjRuH5cuXo6CgAJGRkRgxYgRnsBGpRHm72ADIlWJOoKCaigtFWpfqA9KhQ4fQu3dv+f7UqVMBAKNHj8Zbb72FTZs2AQDCwsIUj9u5cyd69eoFAFi9ejUiIyPRp08faLVaDBs2DEuWLLHNGyCiMpW3iw1gBYmIFSTrUn1A6tWrF4QQJs+Xdk7i4+ODNWvWWLJZRGRBFakgMSBRTceAZF1VfgwSEVV9rCARlR8HaVsXAxIR2R0rSETlxwqSdZkVkI4dO2btdhBRDcYKElH5MSBZl1kBqXPnznjvvffMGu9DRFRerCARlR8DknWZFZDy8/MRHR2Nhx56CBcvXrR2m4iohuE6SETlxzFI1mVWQJo8eTI0Gg3279+P0NBQLF++3NrtIqIahF1sROXHCpJ1mRWQPvroI8TGxqJly5bIzs7GxIkTMWDAAFy/ft3a7SOiGoBdbETlx4UircvsWWzdunXDsWPHMH36dDg4OOD3339HSEgI1q1bZ832EVE1l5+fL/+ir0gFiStpU03FCpJ1lWuav7OzM+bPn48DBw4gJCQEqampePbZZzFy5EjcvXvXWm0kompM6l4DyldBksJUTk6OxdtEVBUwIFlXhVbSfuCBB3D48GF88MEHeOedd7B+/XrExsYiOjq61E+Azz//fIUbSkTVkxSQnJyc4OzsbPbjateurXg8UU3DQdrWVeGtRhwcHDBjxgyMGDECffr0QVJSEiZNmmTyeo1Gw4BERAYqMv5I/3pWkKimYgXJuiq1kvapU6cwduxY/PPPPwCK90Uz9aXT6SzRXiKqZioyxR9gBYmIAcm6KhSQ7t27h9deew3t27fHrl27AAATJkxAZmYmdDqdyS8iopIqMsUfYAWJiAHJusrdxbZz506MGzcOSUlJEEKgWbNmWLFiBR5++GFrtI+IqrmKdrFJFaT8/HwUFhbC0bHCIwaIqiSOQbIusytImZmZGD9+PMLDw+XVtCMjI3H8+HGGIyKqsMoGJIBVJKqZWEGyLrM+cm3ZsgUvv/wykpOTIYRA8+bN8eWXX+Khhx6ydvuIqJqTFnr09PQs1+Nq1aoFjUYDIQRycnLK/Xiiqo4BybrMqiANGTIEycnJAIq3HTl+/DjDERFZhLTQY3kDjkaj4UBtqtG4krZ1md3F1rx5c+zevRsfffQRXF1drdkmIqpBpAqStDJ2eXCgtvVdvXoVd+7cwZUrV9CpUyc899xz9m4S/X8cg2RdZgWkV155BcePH0ePHj2s3R4iqmEqWkECONXfmjIzM3Hu3DkEBgaia9eu6NGjB+Li4rB69WqT3+9bt24hPT3dxi2tucrbxZaSkoKEhARrNqlaMWsM0ocffmjtdhBRDVWZCpIUkFhBsqzMzEw0b94ct27dAgCcO3dOcf7ChQto166dwWN8fX1Ru3ZtBlYb0Q9FhYWFEEJAo9GYvLZhw4YAgCtXriAgIMAmbazKKrVQJBFRZVV0kDbALjZLycnJQbt27fDkk08CABITE+VwZEzJwARArkzk5OTg3r171mkoKZSsGul3uen7/vvv0bhxY/n+2bNnrdqu6oILhxCRXUldbJWpILFiUTm7du3CiRMncOLECdy6dQspKSmlXm8sIOmHort373Ksqg2UDEgFBQVwcnIyuO7ll1/G7du35fv8QGEeVpCIyK5YQbI//XEphw4dkmctm3L+/HmDY/oVp9TUVMs1jkwqWTEyNg7pzp07inAE/PuhhErHgEREdsUKkv0dOnRIvh0XF4fr168DAB555BE89thj8rnWrVsDAA4ePIgJEyZg7ty5EEIAAG7cuCFfd/fuXfn2hQsXcOXKFau2v6YyVkEqSar2BQQE4IknngDAgGQudrERkV1VpoLEQdqWUTIgNWrUCADQo0cPFBYW4o8//gBQHJjOnDmDhIQEuerk4OCA6Oho3Lx5U34OKSDduHED9913H7y8vJCamgqtlp/JLcmcgCSNN2rZsqX8M8aAZB7+ayUiu6pMBYldbOZLS0tDbm6uwfHU1FQkJSXJ9/fu3YvDhw8DABo2bKiY7dSvXz95JpTknXfeQUFBgaKCJHWxScEqPT291EHfVDElA5GxxSIZkCqOAYmI7MoSFSR2sZkmhMD06dPRoEEDNG3aFKtWrcIjjzyCqKgoAMDJkycBAP7+/mjTpg3S0tLkipK/v79cTQKK/8hevHgRV69eRWFhIdzd3ZGXl4dz584ZrSDt2rVLPsZuNsszZwySFJBatGjBgFRODEhEZDf5+fnIy8sDwAqStZw9exYLFixAfn4+UlJSMHr0aOzcuRNLlixBXl4eTp8+DQBo164dNm/eDBcXF/k9JcYFAAAgAElEQVSxDRs2VPx/CQgIQK1atdCoUSM4ODggODgYAHDq1CmDMUhCCGzbtk0+dvXqVWu/1RqnPGOQ9CtI0ocSKh0DEhHZjf4vag7Stg4puHh4eGDChAmKcxcuXMCZM2cAAMHBwWjWrJliUHbDhg3RoUMHODs7o1mzZnIglbRt2xZAcUDSryDdunULBw4cUFSNWEGyPHMCkhRMGzduzApSOTEgEZHdSL+oXV1d4ehY/jkjHKRdNmmKd0hICGJiYvDHH3/IVaLExES5giTNUOvbt6/82AYNGsDT0xMpKSk4ceKEwXO3adMGADBnzhzFOKbly5ejW7duimsZkCzPnICk34XNgFQ+DEhEZDeVGX8EsIvNHFJAqlevHjQaDR599FEMHz4cgDIgSd1l//nPf9CxY0eMHDlSXnSwTp06chjVJwWk0owYMQIA8MEHH+Crr76q/BsiWVkBqaCgQNGFzYBUPgxIRGQ3lZnBBrCLzRz6AUnSqlUrAEB8fDwuXboE4N+A5Obmhri4OKxZs6bM5w4JCSn1/OTJk+XtSwDghRdewIEDB8r3BsiksgZpl+zCZkAqH9UHpNjYWAwePBj+/v7QaDTYuHGj4rwQAm+++SYaNmwIV1dXhIeHGyyDn5qaioiICHh6esLb2xtjx45FVlaWLd8GERlR2QqSu7s7APDn2Yhvv/0WO3bsKDUgbdiwAUDxbDX98+YKCgrCqlWr8OOPP+LSpUvy8wHAe++9h48//hiBgYGKx0yZMgU6nQ53797lH+pKkgKRg4OD4r5E+vlycXGBk5MTA1I5qT4gZWdnIzQ0FMuWLTN6fsGCBViyZAmWL1+OgwcPws3NDX379lWs9xEREYGTJ09i27Zt2LJlC2JjYzF+/HhbvQUiMiE9PR1AxQOSj48PgOLtFOhfR44cwahRo9CnTx85tOgHIKlaJBkzZkyFX2vUqFF48sknERQUJAcvAOjSpQsAKDZJBYADBw7gtddeQ4sWLRASEsJwWwlSIJK6mk0FJKlCy4BUTqIKASB++ukn+b5OpxN+fn7igw8+kI+lpaUJFxcX8d133wkhhDh16pQAIOLi4uRrfvvtN6HRaMS1a9fMfu309HQBQKSnp1vgnRCREELExMQIAGL48OEVenxiYqIAIDw9PS3csqpJp9OJTz75RDz44IMCgOLryy+/VFw3aNAgAUBotVpx+fJli7z+tWvX5NfLzMyUj3/44Yfi008/FZ988olBu+bNm2eR166JXFxcBADh5+cnAIjNmzcrzu/bt08AEE2bNhVCCHHr1i35+15YWGiPJttFRf9+V+mtRpKSkpCSkoLw8HD5mJeXFzp37oz9+/djxIgR2L9/P7y9vdGxY0f5mvDwcGi1Whw8eFDRP64vLy9PHtwGMHETWYPU/VO3bt0KPV6qimRkZCA/Px/Ozs4Wa1tVkZeXh0ceeQQdO3bEU089ZTCVX6JfQdJoNPj+++/x1ltv4b777jPoBqsof39/LF26FD4+PnL3J1DcrQYAOp0OJ0+eVPQILFq0CNOnT+c2JBUgjUGSxuKVXEm7ZAVJf6xfVlYWvLy8bNHMKqtKB6SUlBQAxVNR9TVo0EA+l5KSAl9fX8V5R0dH+Pj4yNcYM3/+fLz99tsWbjER6TM2PqY8vL29odVqodPpcOfOHYNtMGqCrVu3Yt++fdi3bx90Op3inLu7u9yFVfJ77OLigvnz51u8PZGRkSbPabVaxMTEYPLkyXBwcMB9992HO3fu4M6dO6hfv77F21KdCSFQVFQEoHiZDKDsLjYXFxe4uLggLy8PGRkZDEhlYGQ3YebMmUhPT5e/uIYHkeVVNiBptVq5+iQ9V02jv8TBp59+Kt8eO3asYk2jilbprKFly5Zo3ry5HIqSk5Pt3KKqRz8MSRUkUwFJv5onjUOStoMh06p0QPLz8wMAxRL30n3pnJ+fn2KFV6C4LJmamipfY4yLi4u8sJb+AltEZDmVDUj6j62pAUn/91tRURHc3d2Rm5uLL774Ql78Eajc99hapIrf9evX7dySqqc8AUm/a61ly5YAgO3bt1u7iVVelQ5ITZs2hZ+fH/7880/5WEZGBg4ePIiuXbsCALp27Yq0tDR5d2oA2LFjB3Q6HTp37mzzNhPRvxiQyi8/Px+bNm3CqVOnABjucRYWFiavlK2/0ay3t7ftGmkmf39/AMUBadu2bVi0aBHi4uLs3KqqQX8NpPIEpIiICADAqlWrrN3EKk/1Y5CysrJw/vx5+X5SUhKOHj0KHx8fBAUF4ZVXXsG7776LFi1aoGnTppg9ezb8/f3xxBNPACieztqvXz+MGzcOy5cvR0FBASIjIzFixAj5h5OI7IMBybSioiJERUWhVatWKCgoQNu2bdGtWzeEhYXh4sWLuO+++3Du3DmDgHT//ffLt/v37w+gePKKGgdBSxWkY8eOYfz48SgsLISrqysuXbrEMUll0A9D5o5BAoBnnnkGUVFROHbsGM6dO4cWLVrYoLVVk+oD0qFDh9C7d2/5/tSpUwEAo0ePxsqVK/H6668jOzsb48ePR1paGnr06IGtW7eiVq1a8mNWr16NyMhI9OnTB1qtFsOGDcOSJUts/l6I6F9CCIsEpOo6Bumvv/5SzPby8vLC8uXLcfHiRQDA+fPnodPpSg1ITZo0wcmTJ1GnTh3bNLqcpIC0YcMGuSJy7949rFy5EtOmTbNn01RPCkNarVauGJYMSNIAff2A5OPjg+DgYBw/fhwXLlxgQCqF6gNSr169IIQweV6j0WDu3LmYO3euyWt8fHzMWjafiGwnJydHXtC1MgOIq2sF6fjx44r76enpWLRokeLY3bt3DQJS27ZtFffN2S/NXqQqvjRI28nJCQUFBfjss8/w2muvQaPR2LN5qiaFIUdHR3l5C3MqSEDx9/348eMcHF8G9dVciahGkAKNs7OzYpZNeVWXgPTNN9/g77//lu+XDEhAcUVdX0pKisEfuZIBSc1KLsuwdOlSODs748KFC0hKSrJTq6oGKQw5OTnJmwqbG5CksWnXrl2zdjOrNAYkIrKLkrvMV5QUkKrydiN79uzB888/r5g4IgWktWvXYvHixfJxrVYrr/126tQpg8UBq9LYnZLjQB999FGEhYUBgGJT21mzZiE6OrrU3oSaRuqS1A9IZS0UKSlZuSPjGJCIyC4sMf4I+DcQlLbwq9qdOHFCcV+n0yEhIQEA0K5dO/Tr108eZB0ZGSmPG5GmagcGBmLHjh04ffq0DVtdefoVpIceeghNmzaVQ+LBgwcBAFeuXMH8+fMxb968Kvf+rKkyFSQGJPOofgwSEVVPUnm/5Er45RUUFAQAuHz5cqXbZC/SishA8eD1U6dOITs7Gy4uLmjRogUcHR2xd+9eODs7o0OHDhg+fDiA4uoSUFx50Z/MUlUEBgaiX79+0Ol0WLduHTQaDTp37oylS5diz549ACAHRQDYuHGjqsdU2ZL+GCRjAWns2LFyFY4BqWIYkIjILs6cOQMAih3gK0LaLf7u3bvIzMw0+GNQFehvEZKdnY133nkHQHHwcXQs/jUtre0GQN4+SdojcsCAAbZqqkVptVr89ttvimNdunQBABw5cgTPP/88mjdvLp/76aefMGvWLJu2Ua1KqyAJIfDll1/K15oag8SAVDp2sRGRXUjdJcHBwZV6Hg8PD3ka+6VLlyrdLnvQ3xh7//79WL9+PTQaDf773/8avV5/f0lHR0c8+uijVm+jrTRr1gyvvPIKgOKB62+99ZZ87tChQ3K3UU1nbAySFJD0txFp27YtQkJCFI+VKkgpKSmKBSdJiQGJiOzCUgEJ+LeKVFUDUlpamnz7ww8/BAA88cQTaNeundHr9QPSQw89VK22QtJoNFi8eDH27NljdDPVqvr/2NJKqyBJe4fWq1cPCQkJ8kKSEl9fX3mT55JbcdG/GJCIyOZyc3Pladz6+4VVVJMmTQCo94+nTqfD9OnTsW7dOvmYEALr169HfHy84hP/1q1bAQAvvviiyefTD0gDBw60Qovtr0ePHoiKipLvS3/k//nnHwDArl27sHv3bns0TRVKC0jS2liBgYFGH+vg4CDvRcpuNtMYkIjI5s6dOwedTgcvL69SN402l9orSPv27cOCBQswYsQIORj+8MMPeOaZZ9CnTx+kpqYqrq9Xrx769u1r8vn0A1JVHX9kjsmTJ6NevXoICgpCnz59ABQHpOTkZPTu3Ru9evWSFxutaUobpC1VkEwFJIDjkMzBgERENpeYmAiguHpkidWS1R6Q9Fe7fvPNNyGEkMcXpaWl4a+//lJc36RJEzg4OJh8vuDgYLi7u6NDhw4WqcCpVd26dXHy5EnEx8fLSxv8888/+PXXX+VrqvL6V5WhX0EquZK29O8tICDA5OM5k61snMVGRDYnBZlmzZpZ5PmaNm0KoHj/MjXOZNNfsXjNmjUYPHgwjh49avQ8oKwQGVO/fn2cPXsWrq6u1X47Dul7od+NevbsWfn8nTt34OfnBwcHBwgh8PTTT6OwsBA//PCDKjfotZTSBmmbU0GSAhJX0zat+v7rISLVkgKSVPmprPDwcDRu3BhXr17F7NmzLfKclqT/R0in0+GZZ54p9fqyAhJQvMiit7d3pdtWVUj/VuLi4rBt2zb5eGhoKB544AEUFBQgOTkZ33//PTZu3IhTp07Zq6k2YWwMkrSSNitIlsGAREQ2Jy3qKC3yWFnu7u5YuHAhgH9Xly6PgoICzJo1C7///nu5H5ubm4sRI0ZgzJgxisqGPumPkLTAo8RUUDInINU0+hWkkuOOjh8/jt27d+PChQvysbi4OFs2z+bMmcXGMUiVw4BERDYnVZAsFZCAf2fD3bhxo9yPXbJkCebPn49+/fqV+7Fbt27FunXrsHLlSjz99NPycf19w6QK0lNPPYXQ0FD5+OTJkw2mYAMMSMY0bdpUXjQTgMH37bfffsP58+fl+9JWJdWVqUHaQghWkCyEAYmIbE6qIFmqiw34d8uS27dvG+xJZUphYSHeeOMNvPbaa/KxzMxMFBUVYc6cOdixY0eZz3HkyBH59unTp1FUVIRr164hMDBQnqovBaRGjRrhueeek6/v1KmT0XWgKrv9SnXk6emJdevW4fnnn8f777+Pp556SnF+8+bNigpSdQ9IpsYgpaam4t69ewD+rRIZw4BUNg7SJiKbysrKkqe1W7KCVLduXTg4OKCoqAhNmzZFVFQUpk2bVupjNm3aZLBa9YkTJ3Dy5EnMnTsXAIzuIC+EwPLly+Ht7Y1Dhw7Jx/Pz83H16lVER0fj2rVrWLFiBf73v//Jf4QaNWqEyMhInD17Fp06dYKjoyNatmyJ+Ph4xfOzgmTc0KFDMXToUADAlClTFOfOnTuHn376Sb5/4sQJ5OXlwcXFxaZttBVTXWxS91r9+vVRq1Ytk4+XAtLt27er9fepMhiQiMohPz8fsbGx6N69u9GuESqbVD3y8vKy6ArQWq0W9evXR0pKCq5du4bXX38dr776aqkzmYxViI4fP459+/aV+lo7duzAhAkTjJ7btm0bVq9eLd//559/5K1EGjZsCBcXF3z++efy+fvvv9/gORiQyubj42NwTFqdHSjeAPj27dulVlGqMlMBqaxFIiU+Pj5wcXFBXl4erl+/Lo/xon+xi42oHKKjo/Hoo4+WWZkg06zRvSYp2TX1xRdf4IEHHsCiRYuMXi+txNygQQO0bNkSQHFAun79unxNUlIS0tPTFY+LjY01eK7u3bsDAMaNG6c4Li1lUK9ePaOf0idMmIBmzZopliZgQCqbfkCS/t+VpL8AZ35+PrKysqzeLlsxNQZJqiCVNv4IKN7SRbpGreuH2RsDElE5SDOlli1bZueWVF3WGKAtKRmQ/u///g9HjhzBggULFMevX7+Oq1evIiEhAQCQkJAgb4q6a9cuOcQBxQGnYcOG6NmzJ1555RUkJibKFSZpp/m2bduiffv2pbbt+eefN3q8Tp06OH/+PD799FP5WL169cx4tzWbfkB64YUX5AqIt7e3/G9LCkg6nQ5hYWFo27YtMjIybN5WazA1BsncChIAefFNaeFWUmJAIirFgQMHMHz4cMWiflQxBQUF+Prrr7Fr1y4AtglIkps3b8pjic6cOYOAgAD5D0hYWBjq1auHPn36wMvLC6dPnzaYrn/v3j3s2bMHH3/8Me6//3658vTDDz9gw4YNWLdunRyWAODrr79GeHi4fD8iIsJkFQso/jSv0+nk+9LKyGSafkAKCgpCQkICtm/fjgMHDsjdaqmpqZg0aRKeeuopnD59GpcvX8bXX39tryZbVFljkMqqIAFAq1atADAgmcIxSEQmxMfHo2vXrgCKP+VLu6wDxWNJqHyee+45rF+/Xr5vjS62kvu6denSBQcOHAAAjBo1CiNHjkRKSoocRmrVqiVXA319fbFq1So8/vjjBs/r7u6OmJgYrFu3Dr/99huA4llV999/vzxt//jx4/L1gwcPxs6dO+X7ISEhZbb9ySefRIsWLdC7d+/yvOUaSz8g+fn5wc3NTd6vTTp3+vRpxMTEKB73xRdfIDIyssqvQG5sqxFpkgBgXgWJAal0rCARmfDNN9/ItxMTExV/AI3NbKrpMjMzERUVZTDAWQiBPXv2KMIRYJ0Kkv4YnytXrmDv3r3yeJ7Vq1dj0KBBiurQrl270K1bN/n+kCFDMH36dIPnHTlyJEaPHo1ff/0VmzdvRseOHfHqq68q9ksbNGgQunfvjujoaNSpU0fxCb5NmzZltt3d3R2JiYn47LPPyvema6i6devKt0sGYykgGRtsf/z4caSkpFi3cTZQVgWJAanyWEGiGu/mzZtYvHgxoqKiFL9opV3XAeDWrVuKdVVu3boFnU5Xrfd6Kq+5c+diyZIlWLJkiRwgExMT0a1bN3kLBH3WCEiZmZnybSmgBAYG4ubNm/Lxw4cPAwCWLl2Kzp07GzzHu+++i+bNm2PixInyHyH9rrtBgwZh0KBBBo/z8PDA3r175fv6VUZzAhKAKl/VsKU6derIt+vXr684JwWkkpsAS86fP1/lq8CmBmlL4c+c9ycFpKSkJOTn57NrtwT+dqca7+OPP8Z7772HmTNnKo5fvHhRvn3mzBm8/fbb8v2ioiKkpaXZrI1VwYkTJwyOLVq0CKmpqUZnD1mji+0///kPAKBXr17ysZKfpP/8808Apmc+OTo6Yty4cXJ3DVCxWWX6YYdTqC2vTp06mDFjBmbOnGkwqF0KSKZ+RvVX3K6qjA3SzsjIkH/WzPk36+/vD3d3dxQVFSkW2aRiDEhU40m/GH7//Xe58iGEUFSQgOJfPl26dEHt2rUBFFeRaqKEhATF9GljPv30U5w7dw5ubm4mrynZLWIJYWFhSEpKkscJAaa7GkwFJIl+haIiAWn48OHw8fHBM888o+iKI8uZP38+5s2bZ3Dc2BpJ+vQD0pkzZ3DixAnMmDEDUVFRVab7vLTNap2cnMxaY0yj0ciV3OrQ7WhpDEhU40lTuq9fv47jx49j+/btuHbtmvxJTH/H9AkTJsgr0Op329QU27ZtQ0hICBo1aoRVq1YpzulXTCZMmIA2bdoYBAP9X9rWCg1NmjRRrCCsH3QkLi4uZXbxVTYg1a9fH9evX8d3331X7sdS5ZQVkKQPRbm5uejWrRvatWuH999/H0uWLFFUjtXMWECS+Pr6mt1dK1Xfbt++bdkGVgMMSFTjSYMaAeCRRx7Bo48+ilGjRgEo3hpCvwLRvXt3ebxDTawgbd26FUDxH5Yvv/xSca7kL+TCwkKD75G0wGZplSVLM/aHIjg4uMzxY5UNSEDxdH2OK7I9UwFJWqvq/PnzmDNnDrp37467d+8qrjlz5ozi/tWrV1XZJVdaQCo5Jqs0DEimcZA21WiFhYWKzRqlriNprZ6mTZvi3Llz8vmmTZvKv3xqYgVJf98x/dWmAeNBpORGmDNmzICXlxceffRR6zTQiPHjx2PFihUYNmwYbty4ASEEoqKiynwcV7auukwFpM6dO+PIkSM4fPiwPFi/pDNnzmDgwIEAirvag4KCIITAnTt3yqxM2ZI0Bkl/kLakPP9epYBUEz/wlYUBiWq05ORkxQJ9JTVt2hT37t3DjRs3ABSHAOmXT034haLT6bB27Vq0bdsWISEhik1VSwYkab8xffp/hEJCQuDo6IhJkyZZr8FG+Pv74/Lly+Wu5Eh/gADllHJSP/0gExgYKFeJH3zwQXz22WeljjPS388tNTVVvjY+Pl6x+Ke9sYJkfexioxrp4MGDGD9+vDyjqUmTJkZXnu3QoQO++eYb9O7dG3///TeAfz+d6VeWqqs1a9YgIiICYWFhmDFjhmI2WmZmJrKzs+X7xmYMSd0X06ZNw++//279BptQkW4u/YDE5RyqFv2ANHLkSPl2QECAPLVd35EjR+SxYitWrMDEiRMhhMC1a9fka0p2vdlbWWOQzMWAZBp/6qnG2bFjB7p06YL//e9/mDFjBoDiNXkee+wxg2sHDx6M4OBg7NixAw8++CAAoF+/fgCA9evXG8z8qCozYMwldTUCwAcffAAA6NatmzyT79ixY4iIiMCvv/5a6rIHL774YpVbd0aqFhjbYJbUzdvbG507d0ZoaCieeuop+XjdunXxwAMPKK4dM2YM2rVrh9atW8vHPvnkEyQkJCi6iPUXip07dy4iIyNLrT5bGytI1seARDXOmjVr5NvSOKKgoCDMnj0bERERimulndj19ejRA126dEFeXh5mzJghh6Jdu3ahQYMGWLdunRVbbz1nz55FcHAwVq5cKR+Tqmb6wsLC5LAzaNAgrFmzBgMHDix1HZWquPlq9+7dsXv37iozq4n+pdVqceDAAcTHx8v7sgHFA+/1A9KyZcvw5ZdfQqvVGiz7sHfvXqMBKT09HXPmzMGyZcsQFxdn5Xdimv5CkQ4ODooqJwOSZVT5gFRUVITZs2ejadOmcHV1RfPmzfHOO+8oPskLIfDmm2+iYcOGcHV1RXh4eI3oHiHjjJXKg4KC0KRJE3z77bfyZpaLFy82+niNRoN3330XWq0WX3/9tbw1xPDhw3Hr1i2MGDHCeo2vpOTkZPTu3RubN282OLd+/XqcOXNG3mIlKysLJ0+eNLguJCREDkglZwAZ4+DgoFgqoSrp2bOnvKwDVT1arVYRzn19fRUBSb9qVLt2baxevVoeb/bXX38pAtKJEydQVFSkqCTpr64PFM/uNLYoqjXoLxSp/1+AXWyWUuUD0vvvv49PP/0UMTExOH36NN5//30sWLAAS5cula9ZsGABlixZguXLl+PgwYNwc3ND3759kZuba8eWk70Y23dIfyf2UaNG4dKlS6XOdOrTpw/++9//AgBefvllBAYG4s6dO5ZvrIXNmTMHu3btwpAhQwzOHT16FADkAemHDh2CTqdDYGAg7rvvPvm6+++/X9Fdpj/byxhnZ2eO4SG7cXJywunTp5GQkAB3d3d5qj8Ag/FIzz77rDwWae/evYoxSDk5OXjiiSewYsUK+Zi0tUx8fDx27tyJhx56CK1bt0ZGRoY13xIAZRcboBxnxwqShYgqbuDAgeKFF15QHBs6dKiIiIgQQgih0+mEn5+f+OCDD+TzaWlpwsXFRXz33Xdmv056eroAINLT0y3TcLKptLQ0IYQQt2/fFgAMvnbt2lXu50xNTTX6XPb8sbp586bQ6XQmz48ZM0Zu4/z588WVK1fkc82aNRMARN26dYUQQsyYMUMAEE8//bTo2bOn/LjU1FQxadIk+f7ixYtNfh/s/f0gMmbVqlXi888/N3ouIyNDaLVaAUA88MADAoAICQkx+u/az89P6HQ6g+MbN240eN7z58+LI0eOWOw9dOvWTQAQP/74oxBCKF6/PH+nMjIy5MdlZ2dbrH1qUtG/31X+Y123bt3w559/yjt0Hzt2DHv37kX//v0BFG/Cl5KSopie6eXlhc6dO2P//v0mnzcvLw8ZGRmKL6qavvnmG3h7e+Opp56Sp503atRIsZKzfgXJXHXq1FF8GtVXVFRUscZWwvbt2+Hr64tx48aZvEbodT3PnDkTPXr0wO3bt5Geni6Ptblz5w7y8/OxYcMGAMCwYcMUM7rq1KmDnJwc+f6YMWOMvpaxDV2J1GDUqFEmf048PDwQGhoKAPKyFu+88w5eeOEFg2tTUlKM/h0p2TWt0+nQunVrtG/fHpcuXaps8wEYVpCkleEXLVpk1jYjEnd3d3mTWlaRlKp8QJoxYwZGjBiB1q1bw8nJCe3bt8crr7wiD7aVZhnp78Yt3S9t75n58+fDy8tL/jK1nxOpn/QL7Pvvv8dzzz0HoHjsgf4MlIqOM4mJiTE6y8ke+xpt2rQJQPE05e+//97oNVL3meTSpUuYPXu2YlwFAPz666+4cOECXF1dMXDgQDzxxBMA/v0+denSBUDxDC8vLy+jr7Vq1SqMGjXKrtP7iSqiR48eivsBAQFYsWKFYqHUpk2bAgD+97//GTy+5ODtU6dOyR8ytm/fDqB4Usf69esr3Eb9hSIBYMOGDVi7di2mTJlSrufRaDRyl1zJtc1quiofkNavX4/Vq1djzZo1iI+Px9dff42FCxfKA20raubMmUhPT5e/9LejoKpFv/onLe7YunVrRTWlomNkunXrhnv37hkcv3LlCvLy8vDdd9+VubGrpeiHn48++qjMawYPHgyg+NOuNP5IIk2Nfvzxx+Hm5oaoqCh8+umn8piLMWPG4IsvvpArt8Y+QNSpUwerVq0yunwCkZqVDEjSv+8OHTogJiYG8+bNk4OI/qxPiRSQoqKiMHToUMTGxsrn9u7di4yMDPTu3RvPPPMMTp06BaB4Usj48eORm5srV4dKU7KC1KlTJzzzzDMVWvPr/vvvB6BcKZ9Q9QcHBAQEiJiYGMWxd955R7Rq1UoIIcSFCxcEAIO+3549e4rJkyeb/Tocg1R1DRw40GCMwF42aVMAACAASURBVI4dOyw6Rqbk869fv170799fABBTp041uP769eti4sSJIiEhoVyvc/LkSRESEiKPO9AnjZcAIBwcHMTFixfFrFmzxO+//y6eeOIJsX37dtGoUSMBQMTFxYndu3cLAKJFixaKsUnSl6enp7h06ZJZ7XrqqafkxzVu3Fh89NFH5XpfRGpy9epV+d/zmDFjjF5z8eLFUsfdnTx5Ur7dvHlzxbkpU6bIt7/44gtRVFSkOB8QECAyMzNLbWOrVq0EALF79+5Kv9+3335bABDPPvtspZ9LjSr697vKByQfHx/xySefKI7NmzdPtGjRQgjx7yDthQsXyufT09M5SLsG6d69u8EvL51OJ4YMGSIAiIcffrjSr1Hy+SMiIuTb9evXN7j+scceEwBEgwYNyvU6AQEBRkPdzZs3hYeHR5mDpZ2cnAQAcfnyZXH27FkBQLi7u4v27dsLAMLZ2Vm+9vXXXze7XSkpKeKxxx4TP/zwQ7neD5FaTZw4UQwbNkzk5OSYvEb/57HkV3R0dJk/j9JXhw4dDI79/fff4sMPPxQvvfSSKCoqUrxuYmKifN2+ffsq/V63bdsmAIgmTZpU+rnUqMYGpNGjR4tGjRqJLVu2iKSkJPHjjz+KevXqKX65v/fee8Lb21v8/PPP4vjx4+Lxxx8XTZs2Fffu3TP7dRiQqq62bdsKAOLdd98VHTp0ELGxsUKI4lDx/vvvi5SUlEq/Rmkzue6//36D6ytavSr5uNu3bysCoEajES+99FKZv5Bzc3MVs1ekrz59+si3f/3110p/X4iqs759+xr8DEkfMqRZodKXVqsV8+fPN/qBzdjXmjVrhIODgwAg9u/fr3hd/eeOi4ur9PtIT0+XZ+4lJiZW+vnUpsYGpIyMDBEVFSWCgoJErVq1RLNmzUR0dLTIy8uTr9HpdGL27NmiQYMGwsXFRfTp06fc/wgYkKqe/Px8sXfvXuHr6ysAiMOHD1vttYqKisTRo0fFypUrDX7ReXl5GVxfkYB06dIl+TFubm5CCCFmz56teC4XFxeRkJAgWrRoIVeLSn7VqVNHfk43NzdFO0ePHi3f5791otLpd5UBxV1pI0eONPpz16tXL/lxR48eFT///LPBNd27d5cr2yNGjJCPf/HFF/JjU1JSFI+Jj4+3yHvp3bu3AIq73EurmlVFNXaav4eHBz766CNcunQJ9+7dw4ULF/Duu+/K0xaB4lH6c+fORUpKCnJzc7F9+3aDZeWp+omOjkaPHj3k7USsuZqzVqtFaGiovLyEvvT0dMTFxSEqKgrz5s0zWFCysLAQf/75p2LjV2N27twp387OzsZLL72Ejz/+WHFNXl4e2rZti7NnzyIvLw+7du2SZ7lI9Gd0+vn5ybfDwsIU04PLM1WYqCZq06aNfHv16tU4c+YMgoODjV47cOBA+XZoaCiGDBmCAQMGKK5p3LixPFN0y5Yt8vGEhAT59okTJxSPkVb+rqxvvvkGvr6+OHfuHF5++WVMmjTJLsuVqEmVD0hEpkibq0pssd1FySX+pSnwnTp1wpIlSxAdHY1JkyYprlmwYAHCw8Px/PPPl/rc0mwXyWeffSbP0Pvkk0/g7OyMRYsWyec1Gg0efvhh3Lt3T/ELVn8dFv2A9OCDD2LGjBkYOHCgvGQAEZmmH5D8/f3h6Oio2L/R09MTbm5uAIAnn3zS4PG//PILFi5cKN8PCgqSV6nX37JEPxRJS3IEBARgw4YN8vpHldWoUSM8/PDDAICvv/4aMTEx+O233yzy3FUVAxLVGLaqiEgLJ4aHhyMgIMDgfMm1T9544w0AwI8//ljq8yYlJZk899JLLyErKwtTp041OOfo6Ii2bdti1qxZAICxY8fK5/S3DOnSpQv8/f2xZcsWeQkAIjJNPyD5+PgA+Hd9JKB4S6Jjx47hyJEjJhej1d8PLigoyOiabAkJCdixYwdeeuklvPrqqwCAF198EcOHD7fI+5C0a9dOcT8tLU1xPzo6Gh06dKgxCyczIFG1pb8eiIeHh0FXk7XExMTgo48+wurVqxU7iUtrE5UsWwu99ZgkOp0Ohw8fVqxgbWpX+bFjx0Kj0Sg2qzTm3XffRVxcHN577z35mP4+atLij0RkHm9vbzz77LPo2bOnHJb0A1Lv3r3RvHlzhIWFmXwO/f3gAgICFB9apEVob9y4gT59+sgbYwOGYcYSQkJCFPel4QmSefPmIT4+Xt6vrrpjQKJqy9XVVb5tarVna6hduzaioqLg6+urWJjxhRdeKHNByvnz5+PatWtYt24dOnbsiIEDB8orfksBSf8T5ooVKxTdaqXRaDTo2LGjXPIHgKtXr8q39cMcEZln9erV2L17t/wBTL/b+tFHHy3z8U2aNJFv169fX/Hz3blzZ4Mg5OLiAicnJ3Tq1KmSLTdU8rWSk5Pl2/rVJP2FLPPz8xXdgdUJAxJVW/oBqSKry1qCtCJvgwYN4OHhUeaeb7NmzUJoaKi8HcEff/yBmJgYpKWl4e7duwCUY6nGjBlTqfD34osvAjDvFzkRlU2j0SA+Ph67du1SdJ+Z4ujoiJiYGEyePBmdO3dWVJC6deuG0aNHy/cnTpyIf/75B0ePHrXKB5rGjRsr7ksBafny5ejcubN8PDMzU749dOhQNGjQANeuXbN4e+zOOpPqqh9O8696ateubXRquy2lp6eLt956S1y7dk0IIUSbNm3kNj388MMm10Bp3bq1fLtLly4iPj5eABC+vr6iV69eFVpDyZiioiKxb98+kZGRUennIqLKKywslNc/2rRpk7hx44b88/77779b/fUXLFggv17Pnj2FEIYL4UZGRgohhEhISJCPffvtt1ZvW0XV2Gn+RMbk5uYqdpwvawq9tXh6emLOnDly2Vy/gnTfffeZfNyZM2fk24cOHZJnsTRr1gyLFy/Gfffdh2+//bbS7dNqtejatatiLBIR2Y+DgwPCw8PRqFEj9OzZE76+vli1ahXmzJmD8PBwq7/+tGnT5L3jkpOTkZ+fb3CNVFn68ssv5WPm7DlZVFRkdO/KzZs3Y+LEiUZfy55sM2qVyMak7iiJWn7wYmJioNFoMGvWLPz8889lXq/ValFYWChvvtysWTOEhYXh3Llz1m4qEdnJr7/+iqKiInnixahRo2z6+tIHuvPnzxt0uwHFAUkIoZiRW3JD9+vXr+PKlSvyWCkhBIYPH47t27fj0KFDisHpQ4YMAVA8o6/kMij2xAoSVUslP83YawxSSUFBQfj555/RuXPnMscjAcCzzz4LANixYwcAlDobhoiqB61WW+asVGvSHyiekpJicD45ORlJSUmKSR4lA9LIkSPRpUsXHDhwAACwfft2bNy4EVlZWViyZIl8nTQJBQBOnz5tsfdgCQxIVC3pByRfX1/FqrRqMXTo0DKveeSRRxT3O3ToYK3mEBEBUE5wMeb69evYvXu34ljJgHTkyBEIIbBu3ToAUCwvsnbtWuTl5QFQBjBj3W/2xIBE1crJkyfx9NNPY8+ePQCK1/ZJSUkxWNJfDerUqYPAwECD4w888ACA4l9SJbcuad++vU3aRkQ1W0REhOJ+165d8ccff+D/tXfvUVHW+R/A3zPcBblIMoCCsS5peFcKWfKygIrr8Xjb0harpYw0aNVcTatltTTJddXVzGtlnrxsWqi5JREahiIpl0TzBmKgctGQAS9c5/v7gzPPj2dAQ4V5Zpj365w557l8mfMZv+K8/T7f5/sADbf579mzBwAwdOhQAPKApNVqpcUk9+zZAyGEbDXwsrIypKamApCv73a3td6UwoBE7cpLL72EXbt24a233gLQsLqtqVxea05qaiqmTJmCF198EUDDXIOEhASMHTsWSUlJsjVVgIZQRUTU1j777DNs375d2h8yZAhGjBiBzp07A4D0OCJ9kCooKIC/vz/KyspkYenSpUv48ccfce3aNQANSxcADZPBx4wZg59++klqe/78+bb9UPeJAYnajaKiIqSnp8uOmXqg8PX1xY4dO7Bp0yakp6dj48aN8PX1xb59+xASEgLg/x9dwstrRGRMjReO1K+71HhtJ5VKJXvcSW5uLg4ePIiCggLZ+6xZswZAw8r9+lHwrKwsfP3114iNjZXaFRcXm9RjTBiQqN1obvl7Y66g/TDUajWefPJJ2NvbNzn3wQcfID4+Hl988YUClRGRpfL395e2raysAEC21ED//v2lZ9Dpvf/++xg/fjwASKuL6+928/X1lT2KpTmmNIrEgETthv4OiEWLFmHt2rUICAjAM888o3BVD69Dhw544403mr3dloiordja2kr/yRw2bBgAeUDSPylg06ZN0rETJ05IjyKZOHEi1Gq1tO/r6yt7tEpz8vLyWq3+h8V1kKjd0F/j9vDwwPTp0/Hqq68qXBERkXk7c+YMLl++jN69ewMAnnjiCelcv379ADQ8sqhDhw5NJnYPGjQIRUVF0k0z9xpBOnz4MAICAuDu7t4WH+OBcASJ2g19QNJPIiQioofj5eUlC0U2Njb4+OOPERUVhalTp0rHH3vssSY/6+zsLF1uA5oGJHt7e0yZMgVnz57FkCFDTCocAQxIZOaEENI2AxIRUduLiorCxx9/DDs7O+lY4/lKesOHD8e4ceOkfV9fX9nDtqdNm4YdO3bIVtU2JQxIZLYSEhLQuXNn6VZUBiQiImU0viFm1qxZyMnJQc+ePdG9e3cEBQVBpVKhf//+UKlU0Gg0AFq2WK6SVKLxf8HprioqKuDi4gKtVgtnZ2elyyE0LFymX8Y+JycHffr0AQBcv37d5IZqiYjau/j4eBw4cAB79+6VBaaSkhIUFhYiMDAQAHDlyhWcPXsWYWFhRqnrQb+/GZBaiAHJ9Pj4+EjPAoqMjMS2bdtgZWWFmpoaqNUcHCUiogf//ua3CJml69evyx6UqF8g0t3dneGIiIgeGr9JyCxlZmbK9nNzcwFw/hEREbUOBiQyS0eOHAEA9OrVS3acAYmIiFoDAxKZnby8PCxfvhwAEB0dLTvHgERERK2BAYnMzqpVq3D79m0MGzYMsbGxsjvWGJCIiKg1MCCR2UlLSwMAxMTEQK1Wy55Rpl8On4iI6GEwIJFZuXPnDn766ScAwJNPPgkAsrvZGi99T0RE9KAYkMhsJCcno0OHDqirq4NGo4Gvry8AYN68eQCA119/HR07dlSyRCIiaieslS6AqKVWrlwpbffq1QsqlQoA8Nprr2HIkCGyByoSERE9DAYkMhv6tY4AYMKECdK2ra2tdLmNiIioNTAgkVmoqqqSAtL+/fsRERGhcEVERNSeMSCRWThz5gzq6+vRqVMn/OlPf5IurxEREbWFdjFJ+8qVK5g6dSrc3d3h4OCAPn364MSJE9J5IQTi4uLg5eUFBwcHhIeH48KFCwpWTPfr5MmTAIC+ffsyHBERUZsz+4B048YNhISEwMbGBt988w1+/vln/Pvf/4abm5vUZtmyZVi9ejXWr1+P9PR0ODo6YtSoUaiqqlKwcrofjQMSERFRWzP7S2zvv/8+fHx88Mknn0jH/Pz8pG0hBFatWoW3334b48aNAwBs3boVGo0Ge/bswZQpU4xeM90/BiQiIjImsx9B2rdvHwIDA/H000/Dw8MDAwYMwKZNm6Tz+fn5KC4uRnh4uHTMxcUFQUFB0orMzamurkZFRYXsRcphQCIiImMy+4B08eJFrFu3Dv7+/khMTMSMGTPwt7/9DZ9++ikAoLi4GACg0WhkP6fRaKRzzVm6dClcXFykl4+PT9t9CLqnkpISlJaWQqVSoVevXkqXQ0REFsDsA5JOp8PAgQPx3nvvYcCAAYiOjsbLL7+M9evXP9T7LliwAFqtVnoVFha2UsV0v/SjR/7+/ujQoYPC1RARkSUw+4Dk5eWFgIAA2bHHH38cBQUFAABPT08ADaMQjZWUlEjnmmNnZwdnZ2fZi4zvyJEjGDlyJABeXiMiIuMx+4AUEhKCc+fOyY6dP39eesK7n58fPD09kZycLJ2vqKhAeno6goODjVor3b/NmzdL20OHDlWwEiIisiRmH5Bmz56NY8eO4b333kNubi62b9+OjRs3IiYmBgCgUqkwa9YsLF68GPv27UNOTg6ef/55eHt7Y/z48QpXT4Z27tyJLl26ICUlBUDDHDMAiI6OxowZM5QsjYiILIjZ3+b/xBNPICEhAQsWLMA777wDPz8/rFq1CpGRkVKbefPm4datW4iOjkZ5eTmeeuopHDhwAPb29gpWTkIIXL58GUIIfPrpp5g6dSpiY2Px66+/4ssvv8SwYcOQn58PAIiKioK1tdn/dSUiIjOhEkIIpYswBxUVFXBxcYFWq+V8pFaybNkyvPHGG9K+q6srysvLAQBjx47F7t27YW9vDyEEiouLm9yJSERE9Fse9Pvb7C+xkfnat2+fbF8fjoCG9asKCgoghICDgwM8PDyMXR4REVkwBiRSjH7pBFtbW+lY9+7dATQEJP3ltUcffZTPXyMiIqNiQCJF3Lx5U1qK4erVq8jKysKUKVOwd+9eqFQq3Lp1C8ePHwcgf3QMERGRMXDWKxlNXl4eZs+ejbS0NEydOhUA0LlzZ7i7u8Pd3R07duwAAHh7e+PKlStYvnw5AAYkIiIyPo4gkdEsXrwYX331Fa5fv45Vq1YBQJNFPgGgS5cuAIAbN25ArVZj4sSJRq2TiIiIAYlaxbVr1/DRRx+hpqbmrm1+/PHHJscef/zxJscaT8g+dOgQQkNDW6dIIiKiFmJAolbx6quvYtq0aZg+fXqz52/cuIGzZ88CaJhzNHfuXNjb22Ps2LFN2s6ePRuBgYFITEzk6tlERKQIroPUQlwH6d4a32VWX1+Py5cvw9fXFwCwfv16aRVs/fwiAKirq+Pij0RE1Ka4DhIpRgghCzoRERHo1q0bVqxYAQBYsmSJdE7/jDwADEdERGSyGJDooV2/fh11dXXSflJSEgBgzpw5yM7ORkVFhXRu3LhxRq+PiIjofvG/8NQi9fX1UKvVqKmpgZ2dHQBAp9OhoKAAly9fvuvPvfjii1JAWrNmDaKiooxSLxER0cPgCBL9pnXr1sHBwQHjx4+Hvb09Dh06BKDhtn0/Pz8sXLgQADBixAgcOXIEf/nLX7B//34AQFZWFgAgMDAQsbGxcHR0VOQzEBER3Q9O0m4hS56kbfiYj8GDByMtLa3J8djYWKxZs0baHzlypHS57ZVXXsH69evbvlgiIqJGOEmbjKZTp07NHu/Vq5dsf9OmTRg1ahTc3NyklbOJiIjMAecg0W969NFHcenSJWlfq9VCp9M1aRcSEiLb79atGw4cONDW5REREbU6jiDRbyorK5Pt//LLLygpKWnSznAEiYiIyFwxIFETtbW1OH/+PGpra1FdXS3dhRYXFwegYSXsCxcuNPk5tZp/nYiIqH3gNxo1sXjxYvTo0QNdunTBt99+C6BhUce4uDjY2tpCp9Nh2LBhsp+ZP3++EqUSERG1Cc5BoiYOHz4MoOEBtB9++CEA4JFHHoGVlRU0Gg0KCwulthEREYiKisLEiRMVqZWIiKgtMCBRExcvXpS2k5OTAQCdO3cGAFk4AoCePXvimWeeMV5xRERERsBLbCTZunUrNBoNCgoKpGO1tbUAAA8PDwDAvHnzAAChoaGYMGECZs+ebfxCiYiI2hgDkgUQQjS5E+27775DRESE7Pb9F154AaWlpc2+h34EKS4uDj/88AO+++47fPnll/D19W2zuomIiJTCgGQB3nzzTbi7u+PgwYPSsREjRiAxMRGRkZHN/kzv3r1l4Uc/guTo6IinnnqqySraRERE7QkDkgWIj48HAMyZM6fJuaNHjyI9Pb1JUPLw8EDfvn2lff0IEhERkSXgJG0LcrcHxQ4ePLjJsWvXrmHt2rUoKirCtWvXMHbs2LYuj4iIyGQwILVzt2/flrY7duwIoGFOkrW1Nerq6u76c9HR0RgyZAhOnDjR5jUSERGZGgakdi4/P1/a1t+RVlFRcddw5OLigm3btmH06NFGqY+IiMgUMSC1c43XNNLfodbcc9T0evTogTFjxrR5XURERKaMk7RN2Oeff46QkJAmizPq5eTkIC4uDpWVlXd9j7y8PGlbH4yKi4sBAP7+/vjss88QFBQktbGxsWmN0omIiMwaA5IJmzx5Mo4ePYpZs2Y1e/6tt97Cu+++i88//1w6tm3bNvTu3Ru5ubkA5CNI169fR319vRSUPD09ERkZiWPHjklt7hW2iIiILAUDkhm42+KNZ86cAQAsX74cfn5+6NGjB6ZOnYrTp0/j9ddfR21tLf73v/9J7XU6Hfr374+0tDQAgEajafKe+vWOiIiILBnnIJmp2tpaaRXss2fPNjmv1Wrx2Wef4eLFi/Dw8EBlZSXu3LmDU6dO4dSpUwAaRpD09uzZgyVLlmDNmjVGqZ+IiMiUtbsRpPj4eKhUKtllqaqqKsTExMDd3R1OTk6YNGnSPScqmxohRJNjBQUF97xN/+bNm1i8eDGAhuen6XS6Jm0aB6Rx48bhxx9/RM+ePVuhYiIiIvPWrgLS8ePHsWHDBtkK0AAwe/ZsfPXVV9i1axdSUlJw9epVTJw4UaEqW8eFCxfueT4zM1MaPZo+fTqqq6ubtGlugUgiIiJqRwHp5s2biIyMxKZNm+Dm5iYd12q1+Oijj7BixQqEhoZi0KBB+OSTT3D06FHZ5GRTU1VVJW03N4Kkn4Stp1Y335Vz586Fo6Mj/vrXvwJoWE3b1tYWa9asQVhYWOsVTERE1I60m4AUExODMWPGIDw8XHY8IyMDtbW1suM9e/aEr6+vNFm5OdXV1aioqJC9jKm8vFzarq+vBwBs2bIF/fv3x8mTJ5sEJI1Gg9TUVDz//POy49HR0QCAlStXIikpCZWVlbh58yZiY2Pb+BMQERGZr3YxSXvnzp3IzMzE8ePHm5wrLi6Gra0tXF1dZcc1Go20HlBzli5dikWLFrV6rS2l1Wql7fT0dLzyyivYuHEjACA4OBhDhgyRtbe3t0dISAhCQkJw6NAhFBYWYujQoXB2dgYAuLq6SiGRax0RERHdm9mPIBUWFmLmzJnYtm0b7O3tW+19FyxYAK1WK73utlhjW2k8ggRACkdAw/PVEhMTZedtbW2l7V27duHll1/G7t2727ZIIiKidsrsR5AyMjJQWlqKgQMHSsfq6+tx+PBhfPDBB0hMTERNTQ3Ky8tlo0glJSWyu7gM2dnZwc7Ork1rb6ympkYWchqPILVE41GhoKAg2erYREREdH/MfgQpLCwMOTk5yM7Oll6BgYGIjIyUtm1sbJCcnCz9zLlz51BQUIDg4GAFK/9/cXFxcHNzw8GDBwE0jACNGjXqN3/Ozs4Ojz76KADg6aefbssSiYiILIpKNHeLlJkbPnw4+vfvj1WrVgEAZsyYga+//hpbtmyBs7MzXnvtNQDA0aNHW/yeFRUVcHFxgVarleb1tBaVSiVt37lzBw4ODs22c3V1xZIlSxATEwMAePzxx5GYmIjExEQ899xzRh3xIiIiMgcP+v1t9pfYWmLlypVQq9WYNGkSqqurMWrUKHz44YdKlwWg6S38GzZsuGvbGzdu4MaNG1JA0mg08PHxwbRp09q0RiIiIkvTLkeQ2kJbjSAVFRXB29tb2u/atSsuX77cpJ2Tk5P0IFn9iNMf/vAHHDlypNVqISIiam8e9Pvb7OcgmTvDFbGbC0cA4OLiIm3rJ2QPHTq07QojIiKyYBZxic2U6QNS9+7dkZeXd9d2je/Ay8jIwO7duzF37tw2r4+IiMgScQRJYfoVsX/rrrWoqChpu0+fPli0aBGcnJzatDYiIiJLxYCkMP0Ikr+/PyIiIgAAb7/9NpKTk/Hcc8/h9OnT2LVrF2bOnKlkmURERBaFl9gUVlNTAysrK/j7+2PSpEn44YcfMGXKFKjVaoSGhgIAAgICFK6SiIjIsvAuthZqy3WQamtrAfAZaURERK2N6yCZMQYjIiIi08I5SEREREQGGJCIiIiIDDAgERERERlgQCIiIiIywIBEREREZIABiYiIiMgAAxIRERGRAQYkIiIiIgMMSEREREQGGJCIiIiIDDAgERERERlgQCIiIiIywIBEREREZMBa6QLMhRACAFBRUaFwJURERNRS+u9t/fd4SzEgtVBlZSUAwMfHR+FKiIiI6H5VVlbCxcWlxe1V4n4jlYXS6XS4evUqOnbsCJVKpXQ5JqeiogI+Pj4oLCyEs7Oz0uUQ2Cemhv1hWtgfpqUt+0MIgcrKSnh7e0OtbvnMIo4gtZBarUbXrl2VLsPkOTs78x8bE8M+MS3sD9PC/jAtbdUf9zNypMdJ2kREREQGGJCIiIiIDFgtXLhwodJFUPtgZWWF4cOHw9qaV25NBfvEtLA/TAv7w7SYWn9wkjYRERGRAV5iIyIiIjLAgERERERkgAGJiIiIyAADEhEREZEBBiS6p8OHD2Ps2LHw9vaGSqXCnj17ZOeFEIiLi4OXlxccHBwQHh6OCxcuyNqUlZUhMjISzs7OcHV1xUsvvYSbN28a82O0G0uXLsUTTzyBjh07wsPDA+PHj8e5c+dkbaqqqhATEwN3d3c4OTlh0qRJKCkpkbUpKCjAmDFj0KFDB3h4eGDu3Lmoq6sz5kdpF9atW4e+fftKi9sFBwfjm2++kc6zL5QVHx8PlUqFWbNmScfYJ8azcOFCqFQq2atnz57SeVPvCwYkuqdbt26hX79+WLt2bbPnly1bhtWrV2P9+vVIT0+Ho6MjRo0ahaqqKqlNZGQkTp8+jaSkJOzfvx+HDx9GdHS0sT5Cu5KSkoKYmBgcO3YMSUlJqK2txciRI3Hr1i2pzezZs/HVV19h165dSElJwdWrVzFx4kTpfH19PcaMGYOamhocPXoUn376KbZs2YK4uDglPpJZ69q1K+Lj45GRkYETJ04gNDQU48aNw+nTpwGwL5R0/PhxbNiwAX379pUdenleIgAACMVJREFUZ58YV69evVBUVCS9UlNTpXMm3xeCqIUAiISEBGlfp9MJT09P8a9//Us6Vl5eLuzs7MSOHTuEEEL8/PPPAoA4fvy41Oabb74RKpVKXLlyxXjFt1OlpaUCgEhJSRFCNPz529jYiF27dkltzpw5IwCItLQ0IYQQX3/9tVCr1aK4uFhqs27dOuHs7Cyqq6uN+wHaITc3N7F582b2hYIqKyuFv7+/SEpKEsOGDRMzZ84UQvD3w9j++c9/in79+jV7zhz6giNI9MDy8/NRXFyM8PBw6ZiLiwuCgoKQlpYGAEhLS4OrqysCAwOlNuHh4VCr1UhPTzd6ze2NVqsFAHTq1AkAkJGRgdraWlmf9OzZE76+vrI+6dOnDzQajdRm1KhRqKiokEY+6P7V19dj586duHXrFoKDg9kXCoqJicGYMWNkf/YAfz+UcOHCBXh7e+N3v/sdIiMjUVBQAMA8+sI0lqsks1RcXAwAsr+8+n39ueLiYnh4eMjOW1tbo1OnTlIbejA6nQ6zZs1CSEgIevfuDaDhz9vW1haurq6ytoZ90lyf6c/R/cnJyUFwcDCqqqrg5OSEhIQEBAQEIDs7m32hgJ07dyIzMxPHjx9vco6/H8YVFBSELVu2oEePHigqKsKiRYswZMgQnDp1yiz6ggGJyEzFxMTg1KlTsmv6ZHw9evRAdnY2tFotdu/ejRdeeAEpKSlKl2WRCgsLMXPmTCQlJcHe3l7pcize6NGjpe2+ffsiKCgI3bp1w+effw4HBwcFK2sZXmKjB+bp6QkATe46KCkpkc55enqitLRUdr6urg5lZWVSG7p/sbGx2L9/Pw4dOoSuXbtKxz09PVFTU4Py8nJZe8M+aa7P9Ofo/tja2uL3v/89Bg0ahKVLl6Jfv374z3/+w75QQEZGBkpLSzFw4EBYW1vD2toaKSkpWL16NaytraHRaNgnCnJ1dcVjjz2G3Nxcs/j9YECiB+bn5wdPT08kJydLxyoqKpCeno7g4GAAQHBwMMrLy5GRkSG1OXjwIHQ6HYKCgoxes7kTQiA2NhYJCQk4ePAg/Pz8ZOcHDRoEGxsbWZ+cO3cOBQUFsj7JycmRBdekpCQ4OzsjICDAOB+kHdPpdKiurmZfKCAsLAw5OTnIzs6WXoGBgYiMjJS22SfKuXnzJvLy8uDl5WUevx9tPg2czFplZaXIysoSWVlZAoBYsWKFyMrKEr/88osQQoj4+Hjh6uoq9u7dK06ePCnGjRsn/Pz8xJ07d6T3iIiIEAMGDBDp6ekiNTVV+Pv7i2effVapj2TWZsyYIVxcXMT3338vioqKpNft27elNtOnTxe+vr7i4MGD4sSJEyI4OFgEBwdL5+vq6kTv3r3FyJEjRXZ2tjhw4IDo3LmzWLBggRIfyazNnz9fpKSkiPz8fHHy5Ekxf/58oVKpxLfffiuEYF+YgsZ3sQnBPjGmOXPmiO+//17k5+eLI0eOiPDwcPHII4+I0tJSIYTp9wUDEt3ToUOHBIAmrxdeeEEI0XCr/z/+8Q+h0WiEnZ2dCAsLE+fOnZO9x6+//iqeffZZ4eTkJJydnUVUVJSorKxU4NOYv+b6AoD45JNPpDZ37twRr776qnBzcxMdOnQQEyZMEEVFRbL3uXTpkhg9erRwcHAQjzzyiJgzZ46ora018qcxfy+++KLo1q2bsLW1FZ07dxZhYWFSOBKCfWEKDAMS+8R4Jk+eLLy8vIStra3o0qWLmDx5ssjNzZXOm3pfqIQQou3HqYiIiIjMB+cgERERERlgQCIiIiIywIBEREREZIABiYiIiMgAAxIRERGRAQYkIiIiIgMMSEREREQGGJCIyCIIIRAREQGVSgUnJyfk5ubes71Wq0X37t2hUqkwZcoUI1VJRKaCC0USkcUoKipC3759cf36dQQFBSE1NRXW1tbNtv3zn/+ML774At27d0dmZiacnZ2NXC0RKYkjSERkMby8vLB582YAQHp6Ot59991m261ZswZffPEFbG1t8d///pfhiMgCcQSJiCzOK6+8go0bN8LKygqpqakYPHiwdO7EiRMICQlBTU0NVq5ciVmzZilYKREphQGJiCzOrVu3MHDgQJw/fx7du3dHdnY2nJycoNVqMWDAAOTn52Ps2LHYt2+f0qUSkUJ4iY2ILI6joyO2bdsGGxsb5OXl4e9//zsAICYmBvn5+fDx8cGWLVuULZKIFGW1cOHChUoXQURkbN7e3rC2tkZycjIyMzOh1Wqly2779+9Hjx49lC6RiBTES2xEZLF0Oh3++Mc/4vDhw9KxJUuW4M0331SwKiIyBQxIRGTRzpw5g4CAAACAh4cHrly5ctdb/4nIcnAOEhFZtOXLl0vbpaWlsn0islwcQSIii7V9+3ZERkZCpVLhqaeewg8//AAbGxukpaVh0KBBSpdHRAriCBIRWaTc3FxMnz4dQMO6SHv37oWXlxdqa2sRGRmJ27dvK1whESmJAYmILE5NTQ0mT56MyspKPPbYY1i+fDnc3Nzw8ccfQ6VS4dy5c1wgksjCMSARkcWZO3cuMjMzYWNjg+3bt8PR0REAEBERgZiYGADApk2bsGfPHiXLJCIFcQ4SEVmUvXv3Yvz48QCAZcuWYe7cubLzd+7cwaBBg3DmzBm4u7sjJycHXl5eSpRKRApiQCIii1FYWIj+/fujrKwMYWFhSEpKgkqlatIuKysLgwcPRk1NDUaMGIHExMRm2xFR+8VLbERkEerr6/Hss8+irKwM7u7u2Lp1611Dz4ABA/DOO+8AAJKSkrBy5UpjlkpEJoAjSERkEd566y289957AICEhATpMtvd6HQ6hIaGIiUlBXZ2dkhPT0e/fv2MUSoRmQAGJCIiIiIDvMRGREREZIABiYiIiMgAAxIRERGRAQYkIiIiIgMMSEREREQGGJCIiIiIDDAgERERERlgQCIiIiIywIBEREREZIABiYiIiMgAAxIRERGRAQYkIiIiIgMMSEREREQGGJCIiIiIDPwfnc6BTn3ezBUAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x31e10f910>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.text.Text object at 0x3291e7290>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = collect(50:1.0:500)\n",
    "\n",
    "# 5 gaussian peaks\n",
    "p1 = 20.0 .* exp(-log(2) .* ((x-150.0)./15.0).^2)\n",
    "p2 = 100.0 .* exp(-log(2) .* ((x-250.0)./5.0).^2)\n",
    "p3 = 50.0 .* exp(-log(2) .* ((x-450.0)./1.0).^2)\n",
    "p4 = 20.0 .* exp(-log(2) .* ((x-350.0)./30.0).^2)\n",
    "p5 = 30.0 .* exp(-log(2) .* ((x-460.0)./5.0).^2)\n",
    "\n",
    "# some background: a large gaussian + linear \n",
    "bkg = 60.0 .* exp(-log(2) .* ((x-250.0)./200.0).^2) + 0.1.*x\n",
    "\n",
    "# some noise\n",
    "noise = 2.0 * randn!(ones(size(x,1)))\n",
    "\n",
    "# the observation\n",
    "y = p1 + p2 + p3 + p4 + p5 + noise +bkg\n",
    "\n",
    "# making a plot\n",
    "plot(x,y,\"k-\",label=\"signal\")\n",
    "legend()\n",
    "xlabel(\"X\",fontname=\"Arial\",fontsize=18)\n",
    "ylabel(\"Y\",fontname=\"Arial\",fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calling the baseline() function to remove the background\n",
    "\n",
    "see documentation at http://charlesll.github.io/Spectra.jl/stable/PreProcessing/#baseline-subtraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# We define the portions of the spectra where we want to fit the signal background.\n",
    "# roi should be an n x 2 array, see documentation at \n",
    "roi = [0 100.;200 220;280 290;420 430; 480 500]\n",
    "\n",
    "# caling the baseline function with a natural spline that will fit the signal in ROIs\n",
    "y_gcvspl, bas_gcvspl, xy_roi = baseline(x,y,roi,\"gcvspline\",p=1., roi_out = \"yes\")\n",
    "\n",
    "# using the ALS algorithm (Baek et al. 2015), 10^2-10^5 lambda and 0.001-0.1 p values are recommended\n",
    "y_als, bas_als = baseline(x,y,roi,\"als\",p=0.01,lambda=10^6,niter=10)\n",
    "\n",
    "# using the arPLS algorithm (Baek et al. 2015)\n",
    "y_arpls, bas_arpls = baseline(x,y,roi,\"arPLS\",p=0.1,lambda=10.0^6)\n",
    "\n",
    "# using a Whittaker smoother to go through the ROI (Eilers 2003)\n",
    "y_whit, bas_whit = baseline(x,y,roi,\"whittaker\",lambda=10.0^6)\n",
    "\n",
    "figure()\n",
    "\n",
    "# plotting the initial signal and the roi\n",
    "plot(x,y,\"black\",label=\"signal\")\n",
    "scatter(xy_roi[:,1],xy_roi[:,2],color=\"red\",label=\"roi\") # the signal used to fit the baselines\n",
    "\n",
    "# plotting the baselines\n",
    "plot(x,bas_gcvspl,\"cyan\",label=\"gcvspline\")\n",
    "plot(x,bas_als,\"purple\",linestyle=\"--\",label=\"ALS\")\n",
    "plot(x,bas_arpls,\"orange\",linestyle=\"-.\",label=\"arPLS\")\n",
    "plot(x,bas_whit,\"blue\",linestyle=\"-.\",label=\"Whittaker\")\n",
    "legend()\n",
    "\n",
    "xlabel(\"X\",fontname=\"Arial\",fontsize=18)\n",
    "ylabel(\"Y\",fontname=\"Arial\",fontsize=18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.2",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
